I know I billed this as “Ditching ArcGIS” in the poll, but I 1) am not an ArcGIS expert and 2) there may still be things that you must do in ArcGIS to get your work done. If that’s the case, there is the R - ArcGIS Bridge that makes it possible for these two analysis tools to play with one another, BUT I won’t discuss that here…
# Set global knitr chunk options
if (.Platform$OS.type == "unix") {
# Do not specify Cairo device for MacOS
knitr::opts_chunk$set(warning = F, message = F)
} else {
knitr::opts_chunk$set(warning = F, message = F, dev.args = list(type = "cairo"))
}
You may need to install some packages (and their dependencies) that you don’t already have. The chunk below will do that for you automatically using the load_pkgs function.
# List packages required to run the script -------------------------------------
pkgs <- c("tidyverse","ggmap","sf","cowplot","here",
"knitr","ggrepel","rnaturalearth","rgeos",
"FedData","raster","mapview")
# Install and load all CRAN packages provided from a character vector
load_pkgs = function(pkgs) {
new_pkgs = pkgs[!(pkgs %in% installed.packages()[ ,'Package'])]
if (length(new_pkgs) > 0) install.packages(new_pkgs,repos = "http://cran.cnr.berkeley.edu/")
invisible(lapply(pkgs,function(x)
suppressPackageStartupMessages(library(x,character.only = T))))
}
# Load packages
load_pkgs(pkgs)
dir.create(here("Figs"))
dir.create(here("Output"))
Let’s create a very basic “map”, showing the location of acoustic samples along the West Coast. The first uses base graphics (yuck!).
# Load some ship track data
cps <- read_csv(here("Data/cps_nav.csv"))
# A basic
plot(cps$long, cps$lat)
Here are the same data plotted using ggplot2.
# A very basic ggplot2 scatter plot of the same data
ggplot(cps, aes(long, lat)) + geom_point()
Let’s set a few parameters to make the map more spatially accurate. But it’s a long way from a useful “map”…
# Add some more aesthetics and coordinate controls, maybe some labels
# Read locations
locations <- read_csv(here("Data/locations.csv")) %>%
filter(group %in% c("city", "landmark") )
# Refine the "map"
ggplot(cps, aes(long, lat)) +
geom_point() +
# geom_path(aes(group = transect, colour = factor(transect))) +
# geom_point(aes(group = transect, colour = factor(transect), size = cps.nasc)) +
# geom_point(data = locations, aes(long, lat)) +
# geom_text(data = locations, aes(long, lat, label = name)) +
coord_quickmap() +
# coord_map(projection = 'azequalarea') +
# coord_map(projection = 'azequalarea', xlim = c(-128,-123), ylim = c(45,50)) +
theme_bw()
Before we go farther, let’s set some map boundaries based on either our data (cps$lat and cps.long) or manually (mb.lat and mb.long).
# Define lat and long bounds for west coast map
wc.lat <- range(cps$lat) #c(32, 52)
wc.long <- range(cps$long) #c(-130, -116)
# Define lat and long bounds for Monterey Bay
mb.lat <- c(36.5, 37)
mb.long <- c(-122.2, -121.7)
ggmap is a mashup of ggplot and online mapping services. There are numerous map sources available to ggmap (e.g., Google, OpenStreetMaps, Stamen). We’ll look quickly at Stamen and Google maps. Withing each of those, there are multiple options (e.g., Google has satellite, terrain, etc.). Each has different arguments to specify the map boundaries.
Create a basic Stamen map (Toner option). Look, a map (with no data)!
# Set west coast boundaries for stamen maps
wc.bounds.stamen <- c(left = min(wc.long), bottom = min(wc.lat),
right = max(wc.long), top = max(wc.lat))
# Download stamen map of west coast
wc.map.stamen.toner <- get_stamenmap(wc.bounds.stamen, zoom = 6, maptype = "toner-lite") %>%
ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw()
Add our data to the map. Look, a map (with data)!
# Add layers to map
wc.stamen1 <- wc.map.stamen.toner +
geom_point(data = locations, aes(long, lat)) +
geom_text(data = locations, aes(long, lat, label = name)) +
geom_path(data = cps, aes(long, lat, group = transect, colour = factor(transect)),
show.legend = F) +
ggtitle("Basic options")
wc.stamen2 <- wc.map.stamen.toner +
geom_point(data = locations, aes(long, lat)) +
geom_text(data = locations, aes(long, lat, label = name),
colour = "red", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) +
geom_path(data = cps, aes(long, lat, group = transect, colour = factor(transect)),
show.legend = F) +
ggtitle("Formatted labels")
wc.stamen3 <- wc.map.stamen.toner +
geom_point(data = locations, aes(long, lat)) +
geom_text_repel(data = locations, aes(long, lat, label = name),
size = 2, segment.colour = "black", segment.alpha = 0.5) +
geom_path(data = cps, aes(long, lat, group = transect, colour = factor(transect)),
show.legend = F) +
ggtitle("Formatted and repelled lables")
# Combine maps
wc.grid.stamen <- plot_grid(wc.stamen1, wc.stamen2, wc.stamen3, nrow = 1)
# Save map
ggsave(wc.grid.stamen, filename = here("Figs/wc_map_stamen_toner.png"), width = 12, height = 7)
# Print map
include_graphics(here("Figs/wc_map_stamen_toner.png"))
# Set west coast boundaries for stamen maps
wc.bounds.google <- c(mean(wc.long),mean(wc.lat))
# Download stamen map of west coast
wc.map.google <- get_map(location = wc.bounds.google, zoom = 4, source = "google", maptype = "satellite") %>%
ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
xlim(wc.long) + ylim(wc.lat)
wc.google1 <- wc.map.google +
geom_point(data = locations, aes(long, lat), colour = "white") +
geom_text(data = locations, aes(long, lat, label = name),
colour = "white", size = 2, hjust = 0, nudge_x = 0.5, angle = 45) +
ggtitle("Minor formatting")
wc.google2 <- wc.map.google +
geom_point(data = locations, aes(long, lat), colour = "white") +
geom_text_repel(data = locations, aes(long, lat, label = name),
colour = "white", size = 2) +
ggtitle("Repelled labels")
# Arrange plots
wc.grid.google <- plot_grid(wc.google1, wc.google2, nrow = 1)
# Save map
ggsave(wc.grid.google, filename = here("Figs/wc_map_google.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_google.png"))
# Set west coast boundaries for stamen maps
mb.bounds.google <- c(mean(mb.long),mean(mb.lat))
# Download stamen map of west coast
mb.map.google <- get_map(location = mb.bounds.google, zoom = 10,
source = "google", maptype = "satellite") %>%
ggmap() + xlab("Longitude") + ylab("Latitude") + theme_bw() +
xlim(mb.long) + ylim(mb.lat)
mb.google1 <- mb.map.google +
geom_point(data = filter(locations, name != "Monterey Bay"), aes(long, lat), colour = "white") +
geom_text(data = filter(locations, name != "Monterey Bay"), aes(long, lat, label = name),
colour = "white", size = 4, vjust = 0, nudge_y = 0.01)
# Save map
ggsave(filename = here("Figs/mb_map_google.png"), width = 7, height = 7)
# Print map
include_graphics(here("Figs/mb_map_google.png"))
rnaturalearth for shoreline and country files. Below is a comparison of the medium- and large-scale country datasets.
# Import NaturalEarth coastlines and countries
coast_sf <- ne_coastline(scale = "medium", returnclass = "sf")
coast_sf_lg <- ne_coastline(scale = "large", returnclass = "sf")
countries_sf <- ne_countries(scale = "medium", returnclass = "sf")
countries_sf_lg <- ne_countries(scale = "large", returnclass = "sf")
# View unique regions
# sort(unique(countries_sf$region_wb))
# sort(unique(countries_sf$subregion))
# Select only certain regions to speed plotting
na_sf <- filter(countries_sf, subregion %in% c("Northern America","Central America"))
na_sf_lg <- filter(countries_sf_lg, subregion %in% c("Northern America","Central America"))
# Create West Coast map (medium scale)
wc.ne <- ggplot() + geom_sf(data = na_sf) +
geom_point(data = locations, aes(long, lat)) +
geom_text_repel(data = locations, aes(long, lat, label = name),
size = 2, segment.colour = "black", segment.alpha = 0.5) +
scale_x_continuous(name = "Longitude", limits = wc.long) +
scale_y_continuous(name = "Latitude", limits = wc.lat) +
theme_bw() + coord_sf() +
ggtitle("Medium scale")
# Create West Coast map (large scale)
wc.ne.lg <- ggplot() + geom_sf(data = na_sf_lg) +
geom_point(data = locations, aes(long, lat)) +
geom_text_repel(data = locations, aes(long, lat, label = name),
size = 2, segment.colour = "black", segment.alpha = 0.5) +
scale_x_continuous(name = "Longitude", limits = wc.long) +
scale_y_continuous(name = "Latitude", limits = wc.lat) +
theme_bw() + coord_sf() +
ggtitle("Large scale")
# Create map grid
wc.grid.ne <- plot_grid(wc.ne, wc.ne.lg, nrow = 1)
# Save map
ggsave(wc.grid.ne, filename = here("Figs/wc_map_natEarth.png"), width = 8, height = 7)
# Print map
include_graphics(here("Figs/wc_map_natEarth.png"))
mapview(ca_waters, alpha.regions = 0.1) +
mapview(ca_mpas,zcol = "Type") +
mapview(nasc.sf, cex = "bin.level", color = nascPal, alpha = 0.5)
sf packageAccess to the TIGER/U.S. Census line data files
The elevatr package provides easy access to elevation data from a variety of sources. Not fully functional at the moment…
FedData Github repo and tutorial with examples of various data sources (e.g., .
Presently not executed; still under development
# Get a map of Moss Landing
ML.bbox <- get_map("Moss Landing, CA", zoom = 14)
# Plot the map
ggmap(ML.bbox)
# Get the bounding box from its attributes
bb <- attr(ML.bbox, "bb")
# Create an extent polygon from the Monterey Bay bounding box
ML.extent <- polygon_from_extent(raster::extent(bb$ll.lon, bb$ur.lon,
bb$ll.lat, bb$ur.lat),
proj4string="+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
# Map the bounding box to check the results
ggmap(ML.bbox) + geom_polygon(aes(x=long, y=lat), size=3, color='purple',
data=ML.extent, fill=NA)
# Download National Elevation Data
ned_ML <- get_ned(template = ML.extent, label="ned_ML",
res="1", force.redo = F, raw.dir = here("Data"))
# Convert raster to data frame for plotting in ggplot
ned_df <- rasterToPoints(ned_ML) %>%
data.frame() %>%
rename(lat = y,
long = x,
elev = grdn37w122_1)
# Map elevation
ggplot(ned_df, aes(long,lat,fill = elev)) + geom_raster() +
geom_contour(aes(z = elev), alpha = 0.5, binwidth = 1) +
scale_fill_viridis_c(name = "Elevation (m)") +
theme_bw()
Simple Feaures for R: Blogs, presentations, vignettes for the sf package
A Tidy Approach to Spatial Data: Simple Features from our own Eric Anderson
GIS with R: an introduction by Francisco Rodriguez-Sanchez (Pakillo)
An Intro to GIS with R by Jess Sadler (covers sp, sf, and rnaturalearth, among other things).